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Abstract 

Accurate understanding of the subsurface production rate of the radionuclide ^®Ar is necessary for argon dating tech¬ 
niques and noble gas geochemistry of the shallow and the deep Earth, and is also of interest to the WIMP dark matter 
experimental particle physics community. Our new calculations of subsurface production of neutrons, ^’Ne, and ^®Ar 
take advantage of the state-of-the-art reliable tools of nuclear physics to obtain reaction cross sections and spectra 
(TALYS) and to evaluate neutron propagation in rock (MCNP6). We discuss our method and results in relation to pre¬ 
vious studies and show the relative importance of various neutron, ^’Ne, and ^^Ar nucleogenic production channels. 
Uncertainty in nuclear reaction cross sections, which is the major contributor to overall calculation uncertainty, is 
estimated from variability in existing experimental and library data. Depending on selected rock composition, on the 
order of lO^-lO’” a particles are produced in one kilogram of rock per year (order of 1-10^ kg ’ s ’); the number of 
produced neutrons is lower by ~ 6 orders of magnitude, ^’Ne production rate drops by an additional factor of 15-20, 
and another one order of magnitude or more is dropped in production of Ar. Our calculation yields a nucleogenic 
^’Ne/‘’He production ratio of (4.6 + 0.6) x 10“* in Continental Crust and (4.2 + 0.5) x 10“® in Oceanic Crust and 
Depleted Mantle. Calculated Ar production rates span a great range from 29 + 9 atoms kg-rock ’ yr ’ in the K-Th- 
U-enriched Upper Continental Crust to (2.6±0.8)x 10 atoms kg-rock ’ yr ’ in Depleted Upper Mantle. Nucleogenic 
^®Ar production exceeds the cosmogenic production below ~ 700 meters depth and thus, affects radiometric ages of 
groundwater. The ^^Ar chronometer, which hlls in a gap between and ’"’C, is particularly important given the need 
to tap deep reservoirs of ancient drinking water. 
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1. Introduction 

Argon-39 is a noble gas radionuclide with half-life of 269 + 3 years (NNDC, 2016) used in dating of hydrological 
and geological processes with timescales of a few tens to about 1000 years (Loosli, 1983; Ballentine and Burnard, 
2002; Lu et ah, 2014; Yokochi et al., 2012, 2013). Accurate estimates of production rates of noble gas isotopes, 
such as ^®Ar, "^°Ar and ^'Ne, in rocks are necessary for dating groundwaters and in proper interpretation of isotopic 
signatures of gases originating in Earth’s interior (e.g., Graham, 2002; Tucker and Mukhopadhyay, 2014). 

Production of Ar in the Earth’s atmosphere is dominated by the cosmogenic production channel, where Ar is 
produced from the abundant stable "^°Ar (~ 1 wt.% of atmosphere) by the cosmic ray-induced reaction '^°Ar(n, 2n)^® Ar 
(e.g., Loosli and Oeschger, 1968). The notation is the standard shorthand for a transfer reaction ^°Ar H-« —> 2n -H ^®Ar 
and n stands for a neutron. The Ar cosmogenic production keeps the atmospheric Ar/Ar mass ratio at a present-day 
steady-state value of (8.0 + 0.6) x 10“'® (Benetti et al., 2007). 

Below the Earth’s surface ^®Ar can also be produced cosmogenically, by cosmic ray negative muon interactions, 
in particular by negative muon capture on (93.3 % of all potassium; fj." + ^ H- ^®Ar). As the secondary 

cosmic ray muon flux drops off with depth, at depths greater than 2000 meters water equivalent (m.w.e.; ~ 700 m in 
the rock) the nucleogenic production by ^®K(n,p)^®Ar dominates (Mei et al., 2010). 

Argon extracted from several deep natural gas reservoirs shows ^®Ar/Ar ratios below the atmospheric value 
(Acosta-Kane et al., 2008; Xu et al., 2012). Eor this reason, underground gas sources have attracted the attention 
of astroparticle physics projects that require a source of Ar with minimal radioactive Ar for low-energy rare event 
detection of WIMP Dark Matter interactions and neutrinoless double beta decay using liquid argon detectors. 

This paper provides an evaluation of the ^®Ar nucleogenic production rate for specihed rock compositions by 
calculating the rate of ^®Ar atoms produced by naturally occurring fast neutrons. The calculation also yields results 
for other nuclides, in particular the rare stable ^'Ne (0.27 % of neon). We describe in detail the methods of calculation 
of neutron, ^*Ne and ^®Ar yields and production rates, including the use of nuclear physics software tools (SRIM, 
TALYS, MCNP6), and the input data, with the goal to present reproducible results. In section 2 we give the overview of 
the calculation. Section 3 deals with a-particle production. Section 4 discusses production of neutrons, both by (a,n) 
reactions and by spontaneous hssion. Neutron-induced reactions, in particular ^®K(n,p)^®Ar, are described in section 
5. In section 6 we present the results, including coefficients for production rate evaluation (6.3) and an estimate of 
the uncertainty in the calculated results (6.4). A discussion follows in section 7 where we compare our calculations 
with recent results (Yatsevich and Honda, 1997; Mei et al., 2010; Yokochi et al., 2012, 2013, 2014, sec. 7.1), calculate 
the nucleogenic neutron flux out of a rock (7.2), and discuss the geochemical implications (7.3). Concluding remarks 
follow in section 8. 


2. Overview 

Decays of natural radionuclides of U and Th in Earth’s interior produce a particles with specific known energies. 
The a particles slow down (i.e., lose kinetic energy) due to collisions with other atoms that form the silicate rock and 
eventually stop and form “'He atoms by ionizing nearby matter. A small fraction (< 10“®) of the a particles, before 
losing all their energy, participate in (a, n) reactions with light (atomic number Z < 25) target nuclides to produce 
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fast (~ 1 MeV) neutrons and product nuclides. Additional neutrons, about 10% of the total neutrons, are produced 
in spontaneous fission (SF) of (mean SF neutron energy 1.7 MeV). Collectively all of these neutrons can then 
participate in a number of interactions, including scattering and various reactions including (n,p) reactions. We are 
interested in the (n, p) reaction on to produce Ar (Figure 1). 

The number of Ar atoms produced per unit time per unit mass of rock is denoted S 39 Ar (S for source). To 
calculate S 39 Ar it is necessary to know the neutron production rate as well as the neutron energy distribution. The 
neutron spectra are required because reaction cross sections, including that of ^®K(n, p)^®Ar, strongly depend on the 
energy of the incident particles. To calculate (a, n) neutron production and spectra, the production rate of a particles 
and their energy distribution must be known. In sections 3 to 5 we discuss the production of a particles, neutrons, and 
^®Ar. A useful byproduct of the ^®Ar calculation is the production rate of another noble gas nuclide, ^'Ne, which is 
produced in the (a, n) reaction on with a limited contribution from ^"^Mg(n, Q;)^'Ne reactions. 

We evaluate the results for a number of representative rock compositions. Our selection includes the three layers 
of Continental Crust (CC) of Rudnick and Gao (2014), the Bulk Oceanic Crust (OC) of White and Klein (2014), and 
the Depleted Upper Mantle (DM) composition of Salters and Stracke (2004). Table 1 lists the elemental abundances. 
We also provide plug-in formulae allowing inputs of an arbitrary elemental composition to obtain production rates. 


3. Alpha-particle production 


Alpha particles are produced in radioactive decays of naturally occurring ^^^Th, and and their daughter 
nuclides along the decay chains. Other natural a emitters are not included as their contribution is negligible (the next 
most potent a emitter, *^’Sm, gives 20 times fewer a particles per unit mass of rock per unit time than ^^^U). Each 
individual a-decay emits one a particle with specific energies, given by the decay scheme which is characterized by 
the energy levels and intensities. We use decay data—^half-lives, branching ratios, a energy levels, a intensities—from 
“Chart of Nuclides” available at the National Nuclear Data Center website (NNDC, 2016). The resulting a energy 
spectra are shown in Figure 2. Tabulated a energies and intensities as well as charts of decay networks can be found 
in Supplementary Materials. The maximum a energy is 8.78 MeV from decay of ^'^Po in the thorium decay chain. 
The mean a energy in ^^^Th, and chains is 6.0, 6.0, and 5.4 MeV, respectively. 

Decays of ^^^Th, and produce = 6, 7, and 8 a’s per chain. The a production rate Sa —number of a 
particles produced per unit time in unit mass of rock—of each decay chain is proportional to the elemental abundance 
(mass fraction) A of the parent (Th or U) and is obtained by simple multiplication, 

VAl 

Sa = riaAA—^ = riaAC, (1) 

M 

where we have introduced the per-decay-to-rate conversion factor C, 


C = A 


XNa 
M ’ 


( 2 ) 


in which A is the decay constant and the fraction XNaIM is the number of parent nuclide atoms per unit mass of 
element (e.g., number of atoms per 1 kg of natural U) where X is natural isotopic composition (mole fraction), 
M is standard atomic weight, and Na is Avogadro’s number (see Table 2 for values). The three equations (1), one for 
each decay chain, can be recast into one plug-in formula to calculate a production rate from uranium concentration in 
ppm and Th/U ratio. 
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Alpha production rates 5^ are evaluated in Table 3 for the representative rock compositions. 

Each a particle is emitted with a specific initial kinetic energy (Figure 2). It progressively loses its energy, 
mostly by inelastic scattering on atomic electrons and elastic scattering on nuclei. The range of an a particle, i.e., how 
far it travels before stopping, is obtained by integration along its path. 


Range(£'„„) = 


/‘'■"I 


dEg 
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(4) 
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and is a function of the initial energy of the a particle in a given material. The material property -dEdAx is the 
linear stopping power. We use the SRIM software (srim.org; Ziegler et al., 2008), version SRIM-2013.00, to obtain 
the mass stopping power -p^^dEajdx, where p is rock density (listed in Table 1). Mass stopping power is essentially 
identical for all rock compositions we consider and its energy dependence is plotted in Figure 3. The range of an 
average-energy natural a (~ 6 MeV) in the rock is about 25 pm (Figure 3). We assume that a particle propagation is 
isotropic with no channeling in the crystal lattice on the grain scale. 


4. Neutron production 

The overall neutron production rate —number of neutrons produced per unit time in unit mass of rock—in each 

decay chain is calculated from 

S„ = ACT, (5) 

where the neutron yield T, i.e., the number of neutrons produced from decay or hssion of one atom of parent nuclide, 
is the sum of contributions from (a, n) reactions on various target nuclides and from spontaneous fission, 

# targets 

T= 2 T',„ + Tsf. (6) 

i 

We now calculate the individual neutron yields. Energy spectra of the neutrons produced via various production 
pathways are also discussed in the following sections. 


4.1. Neutrons from (a,n) reactions 

Neutron production by (a,n) was comprehensively described and quantified by Feige et al. (1968). Before an a 
particle loses all its kinetic energy, it can enter another atom’s nucleus to form a compound nucleus. The a has to 
have enough energy to overcome the Coulomb barrier, i.e., the electromagnetic repulsion between the target nucleus 
and the a particle. The Coulomb barrier height Vc is estimated from 


T/ _ 1 4142 

47reo r 
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where indices 1 and 2 refer to projectile and target, eo is vacuum permittivity, q is electric change, r is distance, Z is 
atomic number, Ri„, (in femtometers) is the interaction radius taken as the sum of the two atomic radii, which in turn 
are approximated by a common empirical formula assuming nuclear volume to be proportional to the mass number Jl. 
As the Coulomb barrier height is proportional to the charge of the target nucleus, only relatively low Z target nuclides 
allow for compound nucleus formation from natural a particles. Even though the compound nucleus can form with a 
energies below the Coulomb barrier due to quantum tunneling effect, the interaction cross section drops rapidly. 

The compound nucleus is considered a short-lived intermediate state and from this compound nucleus there are 
many possible channels that form various products nuclides. At energies of natural a particles, the most common 
outcome is that the compound nucleus then sheds a neutron to form a product nucleus. This constitutes an {a,n) 
reaction. Eor example, a + form a compound ^^Ne , which then emits a neutron to produce the ^'Ne product or, in 
shorthand, ^®0(a', n)^^Ne. The majority of the relevant (a, n) reactions are endothermic, i.e., their Q value is negative 
{Q value being reactant minus product rest masses), and the incoming a particle has to have energy above a threshold 
Eth for the reaction to proceed. Specifically, the kinetic energy in the center-of-mass system before interaction has to 
be larger than \Q\ if Q < 0. In the laboratory reference frame, this translates into a relation for the threshold energy of 
an endothermic reaction 


Eth - - 


m\ + m 2 


Q, 


( 8 ) 


m2 


where nti is the rest mass of projectile and m 2 that of target nuclide. 

The Coulomb barrier and threshold energy result in a strong energy dependence of the (a, n) reaction cross section. 
Based on the magnitude of cross sections and the abundance of a particular nuclide in natural rocks, we identify 14 
target nuclides that are important for (a, n) neutron production in natural rocks. They are listed in Table 4 together 
with the {a, n) reaction Q values, threshold energies E,h (8), and Coulomb barrier heights Vc (7). 
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The chance that an a particle, emitted with initial energy Ea^, participates in an (a,n) reaction on a light target 
nuclide i is quantified by the thick target—meaning that the a stops within the medium—neutron production function 


PKKo) = N, 


u 
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O-L(^a) 


dEg 

dA' 


dE„ 


(9) 


where N, is the atomic density of nuclide i (# atoms of i per unit volume) and cr'^ „ is the (a, n) cross section for nuclide 
i. To get from neutron production function to (a, n) neutron yield consists of, for each decay chain and target nuclide 
/, accounting for all a decays and a energy levels within each decay, giving 


decays levels 

n,„ = Z 0. Z fkiPiiEul (10) 

A=i ;=i 


where bk accounts for decay chain branching and fki is the a intensity of level / in decay k with a energy Eu. 

To calculate the (a,n) neutron energy spectra, we evaluate the differential neutron production function dPijdEn, 
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( 11 ) 


which is an equation analogous to (9) except that it calls for the neutron spectrum (or differential cross section) 
dcr'^^ldEn, instead of the integrated cross section (£„ denotes neutron energy). To obtain the yield spectrum (or 
differential yield) dY'^^jdEn, we follow the analogy through equation (10). Similar thick target methods for a particles 
have been used to study nuclear reaction rates for astrophysics (Roughton et al., 1983), and these methods have been 
used for thermonuclear reactions in energetic plasmas (Intrator et al., 1981). 

We use the nuclear physics code TALYS (www.talys.eu; Koning et al., 2008), version 1.6 released on December 
23, 2013, to calculate the (a, n) neutron production cross sections cr'^^ as well as the emitted neutron energy spectra 
dcrjj, „ldEn. We use the default TALYS input parameters except for allowing the width fluctuation corrections calculated 
using the Moldauer model at all energies (see the TALYS documentation for details). An example of our TALYS input 
file is provided in Supplementary Materials. The TALYS-calculated (a,n) cross sections for all considered target 
nuclides are plotted in Figure 4. 

The calculated neutron production function P,- (9) of various target nuclides is presented in Figure 5. The figure 
shows the relative importance of various target nuclides, which depends both on the cross section magnitude and the 
nuclide abundance in the rock. Figure 6 shows the neutron yield spectra dT^^/dP,, for each decay chain and target 
nuclide from each of the three actinides. The mean energy of neutrons generated by (a, n) reactions is 1.8MeV in 
^^^Th decay chain, 1.6 MeV in chain, and 1.7 MeV in chain. 


4.2. Neutrons from spontaneous fission 

Spontaneous fission of produces 2.07 neutrons per fission (Shultis and Faw, 2002) and the SF branching ratio 
is 5.5 X 10“^ (NNDC, 2016). The fission neutron spectrum is approximated by the Watt fission spectrum with neu¬ 
tron energy distribution following exp(-£’„/fl) sinh sJbE^ with parameters a - 0.6483 MeV and b = 6.811 MeV 
(SOURCES-4C, 2002), resulting in a mean SF neutron energy of 1.7 MeV. The fission neutron yield spectrum 
dTsF/d£« is included in Figure 6. The number of neutrons produced in spontaneous fission of ^^^Th and as 
well as other nuclides with non-zero SF branching fraction in the decay chains ^^°Th, ^^'Pa), is negligible 

compared to neutrons from (a, n) reactions. 


5. Neutron-induced reactions: {n,p) and (n, a) 

To quantify the last step in the nuclear reaction sequence of Ar production, the ^®K(n, p)^^Ar exothermic reaction 
(Table 4), we use MCNP6, version MCNP6_Beta3, a general-purpose Monte Carlo N-Particle transport code developed 
at Los Alamos National Laboratory (mcnp.lanl.gov). MCNP6 allows us to calculate the ^®Ar yield per one neutron 
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where we specify the neutron energy spectrum for each of the (a, n) and SF neutron production channels. An example 
of our MCNP6 input file is provided in Supplementary Materials. In addition to ^^Ar yields, we also use MCNP6 to 
calculate the ^'Ne yield from the Q')^'Ne endothermic reaction (Table 4). From these nuclide yields and 

neutron production rates discussed in previous section, we can evaluate the ^®Ar and ^*Ne production rates 539Ar and 
S 21Ne- 

6. Results 

6.1. Neutron yields and production rates 

The calculated yields Y and production rates of neutrons produced by both (a, n) reactions and by spontaneous 
fission are reported in Table 5. For the Upper CC compositions we report detailed results for each {a, n) target nuclide 
and each natural decay chain to show the relative importance of various neutron production channels. For the remain¬ 
ing representative compositions, only the '^©(cr, «)^'Ne, SF, and total neutron yields are presented. The maximum 
neutron yields are of the order of 3 x 10“^ neutrons per decay of 1 atom of long-lived radionuclide. With Upper CC 
composition, spontaneous fission is responsible for 11 % of total neutrons produced; with other compositions the SF 
contribution ranges from 8 to 12 %. Of {a, n) target nuclides in Upper CC rocks, ^^Al produces the most neutrons at 
32 % of total, followed by at 23 %, ^^Si at 9.2 %, ^'’Si at 7.8 %, at 7.0 %, ^^Mg at 4.0 % and ^^Mg at 2.4 %. 
These proportions obviously vary with composition, however, these seven target nuclides -H SF neutrons account for 
at least 96 % of neutrons produced for all representative compositions we use. With the magnesium-rich DM com¬ 
position, Mg nuclides alone account for 69 % of neutrons. In terms of contributions of decay chains, with Upper CC 
composition ^^^Th accounts for 57 % and ^^^U for 41 % of neutrons produced. Again, exact proportions vary with 
composition but the trend, ^^^Th > ^^^U » ^^^U, holds. 

6.2. ^^Ar and ^^Ne yields and production rates 

The largest (n, p) and in, a) nuclide yields (per neutron per elemental weight fraction of reacting nuclide, i.e., K or 
Mg) are obtained with highest energy neutron spectra, coming from exothermic (a, n) reactions with positive Q values. 
That is, from {a, n) target nuclides ^^C, '^O, ^®Mg, and in particular ^^Mg, which provides the highest energy neutrons. 
The maximum nuclide yields are of the order of 0.4 per neutron per weight fraction of K for Ar and 0.05 per neutron 
per weight fraction of Mg for ^'Ne. The production rates, per year per kilogram of rock, also factor in the relative 
importance of neutron production channels. Therefore, the seven most neutron producing nuclides (^^Al, ^^Na, ^®Si, 
^“Si, '^O, ^®Mg, ^^Mg) and SF are also responsible for most of ^^Ar production. Neon-21 production by (n, a) is 
dominated by (a, n) neutrons from ^^Mg (55-81 % of ^'Ne produced depending on chosen composition). The results 
for ^®Ar and ^^Ne production rates are reported in Table 5. 

Table 6 summarizes the subsurface production of a particles, neutrons, ^'Ne, and ^®Ar. Depending on selected 
rock composition, on the order of 10^-10'° a particles are produced in one kilogram of rock per year. The number of 
neutrons produced is lower by ~ 6 orders of magnitude, ^^Ne production rate drops by an additional factor of 15-20, 
and production rate of ^^Ar decreases by at least another order of magnitude. In Table 6 we include production rates 
of "^He, ^'Ne, and ^®Ar in units of cm^ STP per year per gram of rock, in order to facilitate easier comparison of our 
results to previous work. 

6.3. Coefficients for plug-in formulce 

In Table 7 we provide coefficients that one can use to calculate subsurface neutron production, ^'Ne production 
by (a, n) and nucleonic Ar production in a rock of arbitrary composition. The word “arbitrary” should be qualified 
with “while close enough to a natural rock composition”. To indicate what “close enough” means, let us state that 
while the provided coefficients are based on a calculation with Upper CC composition, the empirical formulas yield a 
result within 1 % of the actual value for all rock compositions we use here (Upper, Middle, and Lower CC, Bulk OC, 
Depleted Upper Mantle). The one exception is Ar production evaluation for DM composition, where the Upper CC 
coefficients overestimate Ar by 8 %. 

As an example, we use these coefficients to evaluate the neutron and ^®Ar production rate in the {a,n) target 
nuclide channel with a particles from ^^^U decay chain in a rock with Upper CC composition. The appropriate 
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neutron production coefficient in Table 7 is 2.27 x 10* and the neutron production rate is evaluated, 

using elemental abundances from Table 1 , as 

X Ao X Au 

= 2.27 X 10* X 0.480 x 2.7 x 10^'’ (12) 

= 294neutrons/(yearkg-rock), 

which, as a consistency check, agrees with the rate reported in Table 5. With '*0 target, the neutron production rate is 
equal to the ^'Ne production rate by (a ,«). The *^Ar production coefficient ;t39Ai(^^0, ^**U) in Table 7 is 4.68 x 10’ 
and the Ar production rate is evaluated as 

‘539Ar = A'39Ar(*^0, ’**U) X Aq X Ay X Ak 

= 4.68 X 10’ X 0.480 x 2.7 x 10^'’ x 0.0232 (13) 

= 1.41 atoms/(yearkg-rock), 

which again agrees with Table 5 result. A Python script which performs the evaluations (12) and (13) is provided via 
tinyurl. com/ argon39 . 

6.4. Uncertainty in the calculation 

Estimates of rock composition generally carry large uncertainty, especially as one aims to infer the chemistry of 
the deep Earth (for example, the amount of heat-producing elements K, Th, and U in the Earth is only agreed upon 
up to a factor of a few; e.g., McDonough and Sramek, 2014). Still, it is worthwhile to estimate the uncertainty of 
neutron, ’'Ne and *®Ar production, assuming source rock composition is known. The additional uncertainty arises 
from uncertainties in the nuclear physics models we employ and their underlying data sets, i.e., the evaluation of 
stopping power (SRIM), neutron production cross sections (TALYS), and neutron interactions (MCNP6). 

Half-lives of a-decaying nuclides, branching ratios, a energy levels, and a intensities have small uncertainties 
(NNDC, 2016); some uncertainties are listed in Table 2. We estimate the overall uncertainty of calculating a particle 
production to be below 1 %. 

Ziegler et al. (2010) report the accuracy of SRIM stopping power calculations when compared to experimental 
data, which is 3.5 % for stopping of a particles. Since these stopping powers are vital to implementation technologies 
in the semiconductor industry, these values and uncertainty can be used with confidence. 

Estimation of (a, n) and (n, p) cross section uncertainties are challenging. Neutron reactions on common elements 
in the few MeV energy range are important to many applications for fission reactors. Therefore one might expect data, 
models, and codes to have quite small uncertainties. Various nuclear data libraries exist, such as the ENDE/B-VII 
library (www.nndc.bnl.gov/endf) used by MCNP6, which provide nuclear cross section data. The “evaluated nuclear 
data” come from assessment of experimental data combined with nuclear theory modeling. The evaluated library 
datasets come with no uncertainty estimate, however. Often large differences exist between cross sections from 
various libraries and without provided uncertainties, any rigorous comparison is impossible. Various experimental 
datasets are available, some of which include uncertainty estimate. However, some experimental data are at odds with 
cross sections from nuclear data libraries. 

Given this unfortunate situation, we adopt the following strategy for estimating cross section uncertainty. Cross 
sections used in this study are compared to experimental data available via the EXEOR database at www-nds.iaea.org/exfor. 
We integrate each cross section over energy up to a certain upper energy bound and calculate the relative difference 
between the integrated data sets. We use this difference as a measure of the cross section uncertainty (Icr). 

Clearly, this method of uncertainty estimation is not satisfying. However, given the lack of rigorous uncertainty 
estimates provided by existing data libraries, the variability in the cross section data gives us at least some estimate 
of uncertainty. Koning and Rochman (2012) describe the method of rigorous uncertainty estimation in TALYS nuclear 
data evaluations which consists of a Monte Carlo approach to propagation of uncertainty in various nuclear model 
parameters. Unfortunately, the available versions of TENDL (TALYS-based evaluated nuclear data library) do not 
include the uncertainty in {a, n) cross section and angular distribution. 

In the case of (a, ri) reactions, we integrate TALYS-calculated cross sections, used in this study, and the available 
experimental dataset (using linear interpolation between the data points) up to 6 MeV (mean a energy in the ’*’Th 
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chain). The relative difference between the integrals is taken as Icr uncertainty estimate. In the case of the ^®K(n,p) 
reaction, the cross section from ENDF/B-VII.O library used by MCNP6 and an experimental dataset are integrated up 
to 3 MeV (above which the neutron energy spectrum is negligible). Again, the resulting relative difference is adopted 
as Icr uncertainty. 

Table 8 lists the various contributors to the calculation uncertainty. We considered the seven most important 
(a, n) target nuclides, together with spontaneous fission accounting for > 96 % of neutron production. Experimental 
data from Elynn et al. (1978) (EXEOR entires A0509005, A0509007, A0509008) were used for ^^Al, ^^Si and 3°Si; 
Norman et al. (1982) (C0731001) data for ^^Na; combined Bair and Willard (1962) (P0120002) and Hansen et al. 
(1967) (POl 16003) data for '^O. In the case of ^^Mg and ^®Mg where no experimental data are available in the 
EXEOR database, we arbitrarily set the uncertainty at 10%. The experimental ^^K(n,p) cross section dataset is 
constructed from Johnson et al. (1967) and Bass et al. (1964) experimental data. Nolte et al. (2006) report a non-zero 
^®K(n, p)^^Ar cross section in the epithermal energy range. However, the cross section is smaller by at least 4 orders of 
magnitude compared to values for fast (~ 1 MeV) neutrons, and thus will contribute negligibly to the Ar production. 
Plots of cross sections for visual inspection can be found in Supplementary Materials. 

The uncertainty of the neutron propagation calculation is represented solely by the ^'^K{n,p) cross section uncer¬ 
tainty, which is an oversimplification. Strictly speaking, one should consider the uncertainty in cross sections of all 
reactions that the neutron can participate in. Such complete treatment is beyond this study’s scope. However, ^®K(n, p) 
has the largest of all (n, p) cross sections (Khuukhenkhuu et al., 2011), which justifies our simple approach. The N- 
particle simulations performed by MCNP6 also introduce a statistical uncertainty. Argon-39 yield tallies are repeated 
until the statistical uncertainty is below 0.5 %, therefore negligible compared to systematic uncertainty. Statistical 
uncertainty of ^^Ne tallies are as high as several tens of percent in some cases. However, as we show, the contribution 
of ^‘^Mg(n, a)^'Ne production is negligible relative to '*0(a, n)^^Ne in rocks. 

The overall uncertainty of ^®Ar production is calculated using standard error propagation rules following this 
symbolic formula for Ar production rate evaluation. 


[^®Ar production] = [a decay] x [stopping power] * | ^ [(cr, n) reaction] V X [^®K(n,/7) reaction], (14) 

[target j 

where each term in [brackets] carries its own uncertainty and these uncertainties are considered independent. Eval¬ 
uation of neutron and ^'Ne production uncertainty is modified accordingly. We evaluate the overall uncertainty of 
both the neutron production calculation and ^'Ne production calculation at < 20%; the uncertainty estimate of ^®Ar 
production is 30 % (Table 8). We used the composition of Upper Continental Crust in the uncertainty estimation. 
The choice of elemental composition is reflected in the relative importance of various neutron production channels 
which somewhat changes the overall uncertainty with composition. These uncertainty estimates do not include the 
uncertainty in rock composition. To our knowledge, this is the first attempt at consistent uncertainty estimation of the 
^®Ar nucleonic production rate. Yokochi et al. (2014) state a 50% uncertainty of their ^®Ar production calculations, 
however, no reasoning is offered to justify their estimate. 

Our calculation assumes a homogeneous elemental distribution in the rock. This is a good assumption for the 
Earth’s mantle. In the crust, however, the presence of accessory mineral phases introduces a heterogeneity in dis¬ 
tribution of some elements on the grain size scale, typically a few hundred pm and particularly so for K, Th & U. 
In the Continental Crust the distribution of Th & U is heavily controlled by accessory phases (Bea, 1996). Given 
the mean free path of a MeV neutron in the rock, a few centimeters, the geometric effect of this heterogeneity on 
neutron propagation is insignificant. However, as the range of a natural a particle is only few tens of microns, this 
heterogeneity may change the neutron production both in terms of rate and energy content, depending on how various 
(a, n) target nuclides are spatially distributed within the mineral phases (see Figures 5 and 6). This would in turn 
affect the neutron induced reaction yields of ^®Ar and ^'Ne. Martel et al. (1990) developed a simple geometric model 
of spherical accessory phase inclusion in host rock and quantified their effect on neutron yields relative to uniform 
distribution of elements. In the limit of grain size exceeding the most energetic a particle range (25 pm), they find 
that presence of uraninite and monazite inclusions (deficient in the important light element (a, n) targets Al, Na, Si, 
O) in biotite decreases the neutron production by a factor of 5. The neutron and subsequently ^^Ar production clearly 
depends strongly on petrography. 
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7. Discussion 


7.1. Comparison to previous results 

We provide a comparison of nucleogenic production rates of neutrons, ^'Ne, and ^®Ar calculated in this study 
to select recent evaluations (Yatsevich and Honda, 1997; Leya and Wieler, 1999; Mei et al., 2010; Yokochi et al., 
2012, 2014). Similar calculations were previously performed by Fabryka-Martin (1988); Andrews et al. (1989, 1991); 
Lehmann et al. (1993); Lehmann and Purtschert (1997). 

Mei et al. (2010) calculated rates of ^®Ar production for one representative granitic rock composition. We assume 
that their reported Th and U abundance values were erroneously interchanged in their manuscript, as we are not aware 
of any granites where Th to U ratio is < 1. Using their corrected elemental composition, we calculate a neutron 
production rate at 5400 ± 700neutrons/(kgyr), which is identical to the result of 5500neutrons/(kgyr) calculated by 
Mei et al. (2009). The agreement is expected as both Mei et al. (2009) and our calculation use (a, n) cross section 
inputs from TALYS. Mei et al. (2010) then estimate the ^®Ar production rate to be 7 atoms per kg of rock per year. 
Our calculation gives 16 + 5 atoms/(kgyr). The factor of two discrepancy stems from different methods of calculating 
the ^®K(n,p)^®Ar reaction yields. Mei et al. (2010) estimate the ^®Ar production rate as proportional to the ratio of 
^®K(«,p)^®Ar cross section to total neutron absorption cross section. It is not obvious whether—and how— Mei et al. 
(2010) account for the energy spectrum on the neutrons generated by (a, n) reactions and spontaneous fission. In this 
work, we use MCNP6 which calculates large numbers of histories of neutrons distributed along a calculated energy 
spectrum and accounts for all the possible interactions in the material of given composition. Of course, all these 
models assume, perhaps incorrectly, a uniform distribution of elements. 

Yokochi et al. (2012) calculated nucleogenic Ar production for a few representative compositions using a “mod¬ 
ified version of SLOWN2 code”. With their Continental Crust composition with 2.0 % K, 5.0 ppm Th and 1.5 ppm U, 
their calculated production rate is 0.065 ^^Ar atoms/(cm^ yr) which, assuming rock density 2.7 gcm~^, translates to 
24 atoms/(kgyr). When we rescale our result for Middle CC composition which is the closest in K, Th, U abundance 
to their K, Th, U, assuming other elemental abundances unchanged, we arrive at 13 + 4 atoms/(kg yr), another factor 
of two difference in results, in the other direction. 

Yokochi et al. (2013) calculated the ^®Ar production rate for two specific rock compositions. Their rates were 4 to 
6 times above our calculations, in addition to showing some obvious inconsistencies (a lower Ar production rate for 
a composition richer in K, Th, U). Following Sramek et al. (2013) and subsequent discussion, Yokochi et al. (2014) 
updated their result for “Lava Creek Tuff” rock composition to 120 + 60 atoms/(kg yr). Our calculation with their 
K, Th, U input abundances yields 140 + 40 atoms/(kgyr). Their SLOWN2 code calculation assumes mono-energetic 
neutrons at 2 MeV while our calculation accounts for the actual neutron energy spectrum. When we assume all (a, n) 
and SF neutrons are produced with an initial energy of 2 MeV (while keeping the neutron yields unchanged), we 
observe a 35% decrease in the ^®Ar production rate. Yokochi et al. (2014)’s result falls in between our calculation 
with the full neutron spectra (Figure 6) and our calculation with 2 MeV neutrons, and within the large uncertainty 
agrees with our result. 

Table 9 breaks down ^'Ne production into the ^^0{a,n) and ^‘*Mg(n,a') channels. The neutron induced ^*Ne 
production is negligible for crustal compositions and only contributes 3.4% to ^'Ne produced in a Depleted Upper 
Mantle rock. We calculate the nucleogenic ^'Ne/‘*He production ratio at (4.6 + 0.6) x 10“* in Continental Crust and 
at (4.2 ± 0.5) X 10“® in Oceanic Crust and Depleted Mantle. Within uncertainty, our result agrees with Yatsevich and 
Honda (1997) who calculated the ^^Ne/'^He production ratio at 4.5 x 10“®. Their optimistic uncertainty of < 5 % comes 
from the small uncertainty (2%) adopted for the (a, n) yield in their analysis. Using chemical abundances for the 
“crust” and “mantle” compositions of Mason and Moore (1982), which were the input concentrations in both Yatsevich 
and Honda (1997) and Leya and Wieler (1999) studies, we calculate ^'Ne production rates of 469 ±78 atoms/(kgyr) 
((1.77 + 0.29) X 10^2° cm^ STP/(gyr)) for the crust and 4.4 + 0.7 atoms/(kgyr) ((1.67 + 0.28) x crn^ STP/(gyr)) 
for the mantle. Our results lie in between the lower rates of Leya and Wieler (1999) and the higher values of Yatsevich 
and Honda (1997), while both fall within our Icr uncertainty bounds; see also Ballentine and Burnard (2002). 

7.2. Surface flux of nucleogenic and flssiogenic neutrons 

Near the ground, nucleogenic and flssiogenic fast neutrons produced in the rock can be ejected into the atmosphere 
and enhance the near-surface neutron flux derived from cosmic rays. Most of these neutrons that might be ejected 
from the rock would originate within the first few tens of centimeters within the surface and from ~ 3 meter depth 
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at most. To investigate whether the nucleogenic and fissiogenic fast neutrons may be a significant contribution to 
the near-surface neutron flux, we have calculated the neutron flux for two specific rock compositions: granite is 
the median USGS G-2 composition and limestone is the median NIST SRMlc composition (Table 1), both datasets 
obtained from the GeoReM database (georem.mpch-mainz.gwdg.de). The calculated neutron flux out of a flat rock 
surface is 7.8 x 10“^ n cm“^ s“* for granite and 5.1 x 10“® n cm“^ s“* for limestone. The factor of 15 difference reflects 
the disparate Th and U content of the “hot” granite and the “cold” limestone. These fluxes only account for the (a, n) 
and SF neutrons. It turns out that their contribution is minor compared to the cosmic ray neutron production, where 
the flux of 0.4 eV-0.1 MeV neutrons was measured by Yamashita et al. (1966) to be 2.9 x 10“^ n cm“^ s“* at sea level. 
Cosmic ray neutrons therefore constitute the dominant component (> 99 %) of neutron flux at Earth’s surface. 

7.3. Implications for groundwater dating 

Given the consequences of global climate change and the world population increase, a need is arising to explore 
pumping of groundwater from previously unused deep reservoirs. Interestingly, Gleeson et al. (2015) demonstrate that 
only ~ 6 % (1-17 % with uncertainties) of groundwater present at depths down to 2 km around the globe is “modem”, 
i.e., younger than 50 years. This evaluation is based on tritium (^H) dating (half-life 12.32 + 0.02 years). The next 
available radiometric tool on the age scale is ^®Ar (half-life 269 + 3 years), which Alls the gap between the shorter-lived 

and also *^Kr (10.739+ 0.014 years), and the longer-lived (5700 +30 years). The number of ^®Ar dating studies 
and measurements so far were limited due to the large sample size and analytical requirement of the well established 
Low-Level Counting (LLC) method (Loosli and Purtschert, 2005). However, ongoing technological development 
efforts, in particular the Atom Trap Trace Analysis (ATTA; e.g., Lu et al., 2014), now allow for a relatively fast ^®Ar 
age analysis of small sample volumes (Ritterbusch et al., 2014), and are expected to increase the number of analyses 
in the future. Understanding the relationship between the atmosphere- and subsurface-generated ^®Ar is thus essential 
in deep groundwater radiometric dating, in studies of fluid circulation in the Earth’s cmst, circulation pathways and 
residence times of ocean water masses, and in dating of deep ice cores. 

Hydrological and glaciological dating using ^®Ar relies on counting ^®Ar atoms in water/ice samples as ^^Ar ac¬ 
tivity predictably decays exponentially, it is assumed, after sequestration of the sample from the atmosphere. Cosmo- 
genic neutron flux sharply decreases with depth in the soil/rock, and therefore, the atmospheric production mechanism 
^°Ar(n, 2n)^® Ar becomes insignificant. However, any other mechanism of ambient Ar production at depth must be 
factored in the dating analysis in order to obtain correct ages. In principle, several noble gas nuclide production 
mechanisms are possible underground. Cosmogenic production includes neutron spallation, thermal neutron capture, 
negative muon {pT) capture, and fast muon induced reactions with nuclides abundant in the rock (Niedermann, 2002). 
Mei et al. (2010) showed that capture dominates ^®Ar production at depths down to ISOOm.w.e. (~ 700m depth). 
At greater depths, the nucleogenic ^®K(n, p)^® Ar production dominates. 

In most old groundwater, Ar is below detection limit. This is due to the combination of relatively low production 
rate of ^®Ar in rocks, slow release from rock matrix to groundwater, and short half-life. Therefore, ^®Ar studies for 
water resource assessment are typically based on simple radioactive decay where the subsurface production does not 
contribute significantly to the Ar budget and Ar is simply a decay tracer (e.g., Delbart et al., 2014; Edmunds et al., 
2014; Mayer et al., 2014; Visser et al., 2013; Siiltenfufl et al., 2011; Corcho Alvarado et al., 2007). Geochemical 
sites where ^®Ar subsurface production becomes important are high temperature environments such as geothermal 
system (Purtschert et al., 2009; Yokochi et al., 2013, 2014), highly radioactive environments (Andrews et al., 1989), 
or settings with intense water-rock interactions. Yokochi et al. (2012) presented a case study for such high ^®Ar 
settings. 

To account for the nucleogenic ^®Ar in groundwaters, Yokochi et al. (2012) proposed a novel method using mea¬ 
surements of argon isotopic composition in the fluid. The method is based on comparison of a measurement to a 
modeled evolution of ratio in an argon-producing rock ('^'’Ar denotes radiogenic "^°Ar). Its age range ex¬ 

ceeds the usual limit of several half-lives of Ar — after which ^®Ar concentration assumes a steady-state value in 
an Ar-producing rock while "^°Ar continues to accumulate, thus Ar/^°Ar keeps evolving. The input in their model 
are the ^®Ar and "^°Ar production rates in the rock, which then loses argon to fluid. Yokochi et al. (2012) present case 
studies for two sites with available groundwater argon isotopic composition data. Milk River Sandstone and Stripa 
Granite. In our effort to reproduce their results, we have identified a mistake in their groundwater age calculation: 
We are only able to recover Yokochi et al. (2012)’s results if we intentionally divide (instead of multiply) with rock 
density when recalculating from wt.% K to number of atoms per cm^ rock. In Table 10 we show both Yokochi 
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et al. (2012)’s calculation results and our “original recalculation” using their input values and rock density 2.7 g/cm^ 
for both Milk River Sandstone and Stripa Granite, chosen to best match their output. We also show the “corrected 
recalculation,” where we appropriately multiply with rock density when recalculating from wt.% K to number of 
atoms per cm^ rock. Fixing the calculation error results in a downward revision of both the lower and upper limits 
on groundwater age by a factor of ~ 7 (i.e., rock density squared). The corrected upper age limits for Milk River 
Sandstone samples no longer exceed the age of the reservoir lithology. 

For Milk River Sandstone, Yokochi et al. (2012) use a ^®Ar production rate of 0.015 atoms/cm^/yr. Using a 
rock composition of 2.4ppmU, 6.3 ppmTh, 1.07 %K, and 3.91 %Ca (Andrews et al., 1991) we calculate a ^®Ar 
production rate of 9.45 atoms/kg/yr or, using a rock density of 2.6 g/cm^ (Andrews et al., 1991), 0.023 atoms/cm^/yr. 
For Stripa Granite, Yokochi et al. (2012) use ^®Ar production rate of 1.3 atoms/cm^/yr. Using the Stripa Granite 
composition, including 3.84 %K, 33.0 ppm Th, 44.1 ppm U (Andrews et al., 1989), we calculate a ^®Ar production 
rate of 367 atoms/kg/yr or 0.95 atoms/cm^/yr, using a rock density of 2.6g/cm^ (Andrews et al., 1989). Using these 
new Ar production rates, the calculated groundwater ages lie within a factor of 2 of the corrected results of Yokochi 
et al. (2012). The lower limits on groundwater age for the two Milk River Sandstone samples. Well 8 and Well 9, are 
calculated to be 6.1 kyr and 17kyr. The upper limits on age are 22Myr and 27Myr. Both Stripa Granite samples, NI 
and VI, have ^®Ar/"^°Ar above the revised value for the reservoir rock assuming closed system. The upper limits on 
groundwater age of the Stripa Granite samples are calculated to be 0.056 Myr and 0.35 Myr (Table 10). 

8. Conclusions 

Calculations of subsurface nucleogenic ^®Ar production are extremely useful in hydrology applications and guide 
strategies for obtaining argon for particle physics detectors, while ^^Ne production is an integral component of geo¬ 
chemical interpretations of noble gas isotopic composition. We describe and calculate the production of neutrons in 
(a, n) interactions of naturally emitted a particles and the nucleogenic production of noble gas nuclides ^^Ar and ^'Ne. 
We use new nuclear physics codes to evaluate (a, n) reaction cross sections and spectra (TALYS) and to track neutron 
propagation in the rock (MCNP6). 

Nucleogenic ^^Ar production rate in a rock previously calculated by Mei et al. (2010) is a factor of two below 
our result, the calculation of Yokochi et al. (2012) is a factor of two above our calculation, and the Yokochi et al. 
(2014) result lies within our uncertainty margin. We consider the sources of uncertainty in the calculation along 
with identifying areas that require better constraints. The largest uncertainty comes from uncertainties in the nuclear 
reaction cross sections, which we estimate from the variability among available experimental data and cross sections 
from nuclear data libraries. The overall calculation uncertainty is estimated to be ~ 30 % for Ar production, and 
within 20 % for ^'Ne and neutron production. We expect this uncertainty could be decreased if uncertainty estimates 
were available for the (a, ri) cross sections within TALYS and for the ENDF/B-VII library data used by MCNP6. Our 
analysis yields a tighter uncertainty estimate compared to the recent study by Yokochi et al. (2014) who estimated a 
large error bar (50 %) in the calculated ^^Ar production rate. 

The Ar chronometer fills in a gap between and Conventional radiometric ^®Ar dating of hydrological 
reservoirs assumes a simple radioactive decay of initial ^^Ar content in the groundwater. Non-negligible subsurface 
^®Ar production and/or groundwater-rock interaction in specific environments requires a more refined ^®Ar/^°Ar 
chronometer proposed by Yokochi et al. (2012). Their original calculation for the fluid samples from the Milk River 
Sandstone and Stripa Granite reservoirs overestimated the groundwater ages by a uniform factor of 7. An update using 
our new Ar productions rates results in groundwater age limits which are a factor of 4—11 times lower compared to 
the original analysis. 

The Supplementary Materials* document the uncertainty estimate of (a, n) and {n,p) cross sections, and include 
details of natural decay chains, examples of our TALYS and MCNP6 input files, and a link to argon39, a Fortran 90 
code we wrote, which handles the overall calculation. 


* Supplementary Materials are provided as a separate PDF document. 
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Table 1: Composition, as elemental weight fractions (A in g/g), of geochemical units used in this study. Elements with atomic fractions > 10“"^ 
are included (otherwise entered as -) in addition to all of Th, U, and (a, n) target elements, irrespective of their abundances. Continental Crust 
(CC) composition from Rudnick and Gao (2014), Bulk Oceanic Crust (OC) from White and Klein (2014), Depleted Mantle (DM) from Salters 
and Stracke (2004). Granite is the median USGSG-2 composition from GeoReM database. Limestone is the median NISTSRMlc composition 
from GeoReM. Non-zero amount of carbon was included in the CC and OC composition (0.5 %, 0.1 %, and 0.02 % of CO 2 by weight for Upper, 
Middle, and Lower CC, and 0.05 wt% CO 2 for OC). Rock density is taken from CRUST2.0 (Laske et al., 2001) for the crustal layers and from 
PREM (Dziewonski and Anderson, 1981) for the depleted mantle (average from MOHO to 400 km depth). Elemental abundances not provided in 
a particular compositional estimate entered as n/a. 


z 

symbol 

Upper CC 

Middle CC 

Lower CC 

Bulk OC 

DM 

Granite 

Limestone 

3 

Li 

- 

- 

- 

- 

- 

3.40E-05 

n/a 

6 

C 

1.36E-03 

2.73E-04 

5.46E-05 

1.36E-04 

L37E-05 

2.30E-04 

0.110 

7 

N 

8.30E-05 

n/a 

- 

n/a 

- 

- 

n/a 

8 

0 

0.480 

0.469 

0.452 

0.446 

0.440 

0.494 

0.488 

9 

F 

5.57E-04 

5.24E-04 

5.70E-04 

n/a 

LlOE-05 

1.28E-03 

n/a 

11 

Na 

0.0243 

0.0251 

0.0197 

0.0164 

2.15E-03 

0.0300 

L88E-04 

12 

Mg 

0.0150 

0.0216 

0.0437 

0.0621 

0.230 

4.72E-03 

2.50E-03 

13 

Al 

0.0815 

0.0794 

0.0894 

0.0831 

0.0227 

0.0815 

3.40E-03 

14 

Si 

0.311 

0.297 

0.250 

0.234 

0.210 

0.318 

0.0320 

15 

P 

6.55E-04 

6.55E-04 

4.36E-04 

4.36E-04 

- 

6.11E-04 

L75E-04 

16 

S 

- 

- 

3.45E-04 

n/a 

- 

- 

n/a 

17 

Cl 

3.70E-04 

1.82E-04 

2.50E-04 

n/a 

- 

- 

n/a 

19 

K 

0.0232 

0.0191 

5.06E-03 

6.51E-04 

6.00E-05 

0.0370 

2.40E-03 

20 

Ca 

0.0257 

0.0375 

0.0685 

0.0843 

0.0250 

0.0130 

0.359 

22 

Ti 

3.84E-03 

4.14E-03 

4.91E-03 

6.59E-03 

7.98E-04 

2.84E-03 

4.20E-04 

24 

Cl- 

- 

- 


3.17E-04 

2.50E-03 

- 

- 

25 

Mn 

7.74E-04 

7.74E-04 

7.74E-04 

8.52E-04 

L05E-03 

- 

- 

26 

Fe 

0.0385 

0.0460 

0.0655 

0.0635 

0.0627 

0.0185 

2.88E-03 

28 

Ni 

- 


- 


L96E-03 

- 

_ 

38 

Sr 

- 

- 

_ 


- 

4.75E-04 


56 

Ba 

_ 

- 

- 


_ 

L88E-03 


90 

Th 

1.05E-05 

6.5E-06 

L2E-06 

2.10E-07 

L37E-08 

2.47E-05 

9.49E-07 

92 

U 

2.7E-06 

1.3E-06 

2E-07 

7.00E-08 

4.70E-09 

L95E-06 

L42E-06 

K/U 


8609 

14687 

25320 

9300 

12766 

18974 

1690 

Th/U 


3.9 

5.0 

6.0 

3.0 

2.9 

13 

0.67 

p in gem ^ 

2.70 

2.88 

3.05 

2.83 

3.42 

2.7 

2.5 


Table 2: Atomic and decay quantities for natural long-lived parent nuclides. Natural isotopic composition X and relative atomic mass M from NIST 
(2016), half-lives rj /2 from NNDC (2016), A = ln2/ti/2. Atomic mass truncated to six digits. Numbers in parentheses give uncertainty in last 
digits, otherwise uncertainty beyond shown digits. The conversion factor C (eqn. 2) is the number of decays per unit time per unit mass of parent 
element._ 


Quantity 

Symbol 

Unit 

2J2Th 

-T35u 


Elemental abundance 

A 

kg-elem/kg-rock 




Natural isotopic composition 

Z 

mol-nucl/mol-elem 

1.0000 

0.007204(6) 

0.992742(10) 

Standard atomic weight 

M 

gmoL' 

232.038 

238.029 

238.029 

Half-life 

hji 

10 '^ yr 

14.0(1) 

0.704(1) 

4.468(3) 

Decay constant 

A 

io-'^-‘ 

1.57(1) 

31.20(4) 

4.916(3) 

Number of a’s per decay chain 

na 


6 

7 

8 

Per-decay-to-rate conversion 

C 

10 ® kg-elem“' s“' 

4.072(29) 

0.5759(9) 

12.35(1) 
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Table 3: Calculated production rates Sa of o; particles, as number of or’s produced per second in 1 kilogram of rock, in a particular decay chain 
(columns 2-4) and the total (column 5). 


Composition 

2 j2Th 

225U 


Total 

Upper Continental Crust 

256 

10.9 

267 

533 

Middle Continental Crust 

158 

5.24 

128 

292 

Lower Continental Crust 

29.2 

0.806 

19.8 

49.8 

Bulk Oceanic Crust 

5.11 

0.282 

6.91 

12.3 

Depleted Mantle 

0.334 

0.0190 

0.464 

0.817 


Table 4: Tai'get and product nuclides, Q values, threshold energies E^h (eqn. 8), and Coulomb barriers Vq (eqn. 7) of (n, n) reactions. Also Q and 
Efh of neutron induced reactions. In cases where the product nuclide is unstable*, we show the final stable nuclide as well. Energies in MeV. 
Nuclide rest masses for calculation of Q and £,/, taken from NIST (2016). 


Target 

Reaction 

Product 

<2 

Eth 

Vc 

27 A1 

(a, n) 

JOp- ^ 30sj 

-2.6425 

3.0345 

6.8012 

23Na 

(a, n) 

26 Al* ^ 26]y[g 

-2.9659 

3.4823 

5.9577 

29Si 

(a, n) 

32s 

-1.5258 

1.7365 

7.2107 

20Si 

(a, n) 

33s 

-3.4933 

3.9598 

7.1571 

180 

{a, n) 

2'Ne 

-0.6961 

0.851 

4.5626 

2®Mg 

{a, n) 

2‘>Si 

0.0341 

_ 

6.3298 

25Mg 

(a, n) 

23Si 

2.6536 

- 

6.3838 

19p 

(a, n) 

22Na* ^ 22Ne 

-1.9513 

2.3624 

5.0754 

>70 

{a, n) 

20Ne 

0.5867 

- 

4.6168 

sepe 

(a, n) 

59Ni* ^ 59 co 

-5.0961 

5.4607 

11.5272 

41k 

(a, n) 

44Sc* ^ 44ca 

-3.3894 

3.7206 

9.0555 

48Ti 

(a, n) 

5icr* ^ 5iv 

-2.6853 

2.9094 

10.1118 

'^c 

(a, n) 

16 o 

2.2156 

_ 

3.656 

'‘^Ca 

(a, n) 

47Ti 

-2.1825 

2.3812 

9.3791 


(n,p) 


0.2176 

- 

n/a 

24Mg 

(n, a) 

21^6 

-2.5554 

2.6628 

n/a 
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Table 6: Summary of production rates of "^He, neutrons, ^^Ne, and Ar in units of both the number of atoms or neutrons per year per kilogram of 
rock and cm^ STP per year per gram of rock. Compositional estimates are described in Table 1. 


Composition 

# of atoms (neutrons) yr ^ 
'^He neutrons ^'Ne 

kg ^ 

3®Ar 

^He 

cm"* STP yr ^ g 

1 

2®Ar 

Upper Continental Crust 

1.64 X lO^'^ 

10680 

753 

28.7 

6.20 X 10-^^ 

2.84 X 10-^® 

1.08 X 10-2' 

Middle Continental Crust 

8.98 X 10'* 

6114 

416 

13.9 

3.39X 10-'3 

1.57 X 10-20 

5.25 X 10-22 

Lower Continental Crust 

1.53 X 10® 

1129 

70.2 

0.749 

5.79 X 10-'"^ 

2.65 X 10-2' 

2.82 X 10-22 

Bulk Continental Crust 

9.43 X 10® 

6253 

433 

15.3 

3.56 X 10-'^ 

1.63 X 10-20 

5.75 X 10-22 

Bulk Oceanic Crust 

3.79 X 10*^ 

260 

15.8 

0.0235 

1.43 X 10-'4 

5.96 X 10-22 

8.87 X 10-22 

Depleted Upper Mantle 

2.51 X 10^ 

22.4 

1.06 

0.000257 

9.49 X 10-'® 

4.01 X 10-22 

9.68 X 10-22 


Table 7: Coefficients for evaluation of neutron production (columns 3-5) and ^^Ar production (columns 6-8) for arbitrary composition. Second 
column and third row indicate the elements whose abundances (A as weight fractions) cross-multiply a particular coefficient. Neutron production 
coefficients in units of “neutrons per year per kg-rock per wt-frac-target-elem per wt-frac-chain-parent-elem”. Argon-39 production coefficients in 
units of Ar atoms per year per kg-rock per wt-frac-target-elem per wt-frac-chain-parent-elem per K-wt-frac”. See section 6.3 for examples. 




Neutron production 

^^Ar production 


Chain 

22 ^Th 

ZJ5U 

2i8u 

22 ^Th 

2"U 

zjsu 

{a, ri) target 

A 

Th 

u 

U 

Th 

u 

U 

27^1 

A1 

2.65e+9 

3.31e+8 

5.03e+9 

2.59e+8 

2.36e+7 

3.92e+8 

22 Na 

Na 

6.07e+9 

8 .02e+8 

1.23e+10 

5.04e+8 

4.46e+7 

8 .02e+8 

2®Si 

Si 

1.95e+8 

2.53e+7 

3.91e+8 

3.19e+7 

4.05e+6 

6.46e+7 

20Si 

Si 

1 .68e+8 

2.05e+7 

3.17e+8 

1.37e+7 

1.35e+6 

2.51e+7 

'Sq 

O 

8.77e+7 

1.33e+7 

2.27e+8 

1.97e+7 

2.89e+6 

4.68e+7 

2®Mg 

Mg 

1.72e+9 

2.41e+8 

3.72e+9 

3.92e+8 

5.43e+7 

8.13e+8 

25Mg 

Mg 

l.Ole+9 

1.44e+8 

2.22e+9 

2.75e+8 

4.04e+7 

6.33e+8 

19p 

F 

1.60e+10 

2.34e+9 

3.75e+10 

1.78e+9 

2.06e+8 

3.17e+9 

'20 

O 

9.50e+6 

1.43e+6 

2.47e+7 

2.34e+6 

3.39e+5 

5.70e+6 

sepe 

Fe 

1.28e+8 

3.35e+6 

9.55e+7 

6.81e+6 

2.53e+4 

1.62e+6 

41k 

K 

l.lOe+8 

8.92e+6 

1.64e+8 

1.06e+7 

4.72e+5 

1.22e+7 

48Ti 

Ti 

4.35e+8 

2.20e+7 

5.00e+8 

5.88e+7 

2.29e+6 

5.89e+7 

' 2 c 

C 

3.61e+8 

5.48e+7 

9.96e+8 

1.31e+8 

2.10e+7 

4.02e+8 

44Ca 

Ca 

2.98e+7 

2 .22e+6 

4.29e+7 

3.96e+6 

2.13e+5 

4.80e+6 

SF 

1 

3.01e+3 

2.37e+3 

4.44e+8 

2 .88e+2 

3.09e+2 

4.73e+7 
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Table 8: Estimate of calculation uncertainty and its various contributions, assuming chemical composition is known precisely. Contribution of 
the specific (o', n) channel and spontaneous fission to neutron and ^^Ar production is shown in columns 3 and 4. Calculated with Upper CC 
composition. ^ No experimental data available, uncertainty was arbitrarily set at 10 %. 



Uncert. est. 

Neutron 



% 

% contrib. 

% contrib. 

Decay data, a production 

< 1 



Stopping power 

3.5 



^^Al(a, n) cross section 

36 

32 

25 

^^Na(a, n) cross section 

7.7 

23 

15 

®Si(a, n) cross section 

7.3 

9.2 

13 

^°Si(a, n) cross section 

20 

7.8 

5.4 

'^0 (q', n) cross section 

17 

7.0 

13 

^®Mg(a, n) cross section 

lO'f 

4.0 

7.8 

n) cross section 

lO'f 

2.4 

5.7 

Spontaneous fission 

1 

11 

10 

Overall (a, n), neutron production 

12 



Overall {a, n), Ar production 

10 



^®K(n, p) cross section 

28 



Neutron production calculation 

13 



^'Ne production calculation 

17 



Ar production calculation 

30 




Table 9: Calculated production rates of ^^Ne by ia,n) and (n, a) and ^'Ne/'^He ratio. Compositional estimates are described in Table 1. 


Composition 

prod, in atoms/kg-yr 
(a, n) {n, a) Total 

% contrib. 

(a, n) (n, a) 

^'Ne/He 

Upper Continental Crust 

753 

0.159 

753 

99.98 

0.02 

4.59 X 10-*^ 

Middle Continental Crust 

415 

0.165 

416 

99.96 

0.04 

4.63 X 10-*^ 

Lower Continental Crust 

70.1 

0.104 

70.2 

99.85 

0.15 

4.58 X 10-*^ 

Bulk Oceanic Crust 

15.8 

0.0452 

15.8 

99.71 

0.29 

4.17 X 10-*^ 

Depleted Upper Mantle 

1.03 

0.0366 

1.06 

96.55 

3.45 

4.23 X 10-*^ 
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Table 10: Update of Yokochi et al. (2012)’s Table 2, using their notation. Original recalculation: We are only able to recover Yokochi et al. 
(2012)’s results if we intentionally divide (instead of multiply) with rock density when recalculating from wt.% K to number of atoms per 
cm^ rock. We use Yokochi et al. (2012)’s input values and rock density 2.7 g/cm^ for both Milk River Sandstone and Stripa Granite. Corrected 
recalculation: Appropriately, we multiply with rock density when recalculating from wt.% K to number of atoms per cm^ rock. Recalculation 
with new input: We use our new values of ^^Ar production rate and rock density of 2.6 g/cm^ for both Milk River Sandstone (Andrews et al., 
1991) and Stripa Granite (Andrews et al., 1989). The ratio Pr (defined as Pr = ^^Ar/'^^Ar* in the rock; Pp is the ratio in the fluid) is expressed in 
units of atmospheric P^ = 8.1 x 10“^^; a is the argon loss constant; 0 r is porosity of the rock; Agejow and AgCup are the lower and upper limits on 
the groudwater age. 




INPUT 



CALCULATION 



Rock age 
Myr 

K 

wt.% 

39 p 

at./cm^/yr 

Rr closed 

Ra 

-^S^low 

kyr 

al<pK 

year“* 

Age^ip 

Myr 

Milk River Sandstone Well 8 

Yokochi et al. (2012) values 

83.5 

1.1 

0.015 

193 

27.1 

7.3e-4 

118 

Original recalculation 

” 

” 


192 

25.5 

7.2e-4 

107 

Corrected recalculation 


” 

” 

26.4 

3.5 

7.2e-4 

14.7 

Recalculation with new input 


” 

0.023 

45.7 

6.1 

4.3e-4 

21.9 

Milk River Sandstone Well 9 

Yokochi et al. (2012) values 


1.1 

0.015 

193 

75.2 

3.2e-4 

328 

Original recalculation 

” 

” 


192 

71.3 

3.1e-4 

299 

Corrected recalculation 


” 

” 

26.4 

9.8 

3.1e-4 

41.0 

Recalculation with new input 

” 

” 

0.023 

45.7 

16.9 

L9e-4 

26.7 

Stripa Granite NI 

Yokochi et al. (2012) values 

1650 

3.8 

1.3 

155 

Rp > Rk 

L3e-3 

0.45 

Original recalculation 

” 

” 


153 

Rp > Rr 

L3e-3 

0.39 

Corrected recalculation 

” 

” 

” 

20.9 

Rp > Rr 

L3e-3 

0.054 

Recalculation with new input 


” 

0.95 

15.9 

Rp > Rr 

L7e-3 

0.056 

Stripa Granite VI 

Yokochi et al. (2012) values 


3.8 

1.3 

155 

3.3 

2.9e-3 

2.5 

Original recalculation 

” 

” 


153 

3.1 

2.9e-4 

2.5 

Corrected recalculation 


” 


20.9 

0.43 

2.9e-4 

0.34 

Recalculation with new input 


” 

0.95 

15.9 

Rp > Rr 

4.0e-4 

0.35 
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Figure 1: Overview of nucleogenic ^^Ar production. 
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Figure 2: Natural a emission energy spectra. Blue spikes represent intensity of each alpha emitted. Spike intensities sum up to number of alphas 
emitted in particular decay chain (6, 7, 8 for ^^^Th, respectively). Red dashed line shows the mean energy of emitted a particles. 
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Figure 3: Mass stopping power -p~^dEldx for an a particle in rock (dashed line, left vertical axis) and range of a particle (eqn. 4, solid lines, 
right vertical axis) as a function of a energy. Range plotted up to maximum energy of a natural a particle (8.78 MeV from decay of ^^^Po). The 
dilference between range in Upper Cont. Crust (red) and Depleted Mantle (green) simply reflects the difference in rock density (Table 1). Curves 
for Middle and Lower CC and OC fall between the two shown. 
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Figure 4: Energy dependent (a,n) cross sections for target nuclides considered in this study as calculated by TALYS version 1.6. One mbam 

= 10^27 
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Figure 5: Thick target (rr, n) neutron production function P (i.e., neutron yield per one a pailicle; eqn. 9) plotted against the initial a energy for 
various target nuclides. Curves for the six most neutron producing targets are labeled and color-coded. Remaining targets plotted as thin black lines. 
The Upper CC rock composition is used. Different amplitudes for various target nuclides stem from the combined effect of elemental abundance, 
natural isotopic composition, and energy dependence of (cr, n) cross section. 
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yield spectra dYa,nl^E,,i plotted against neutron energy for each decay chain and six most neutron producing targets. The 
3ctrum dFsp/d-En is also plotted (non-negligible only in decay). The Upper CC rock composition is used. Area below 
the neutron yield Ya,n or f'sF (i-O-. neutrons produced per decay of one atom of parent nuclide). 
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